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We present numerical evidence that in a system of interacting bosons there exists a correspon- 
dence between the spectral properties of the exact quantum Hamiltonian and the dynamical chaos of 
the associated mean field evolution. This correspondence, analogous to the usual quantum-classical 
\^ \ correspondence, is related to the formal parallel between the second quantization of the mean field, 

. which generates the exact dynamics of the quantum TV-body system, and the first quantization of 

^\ ' classical canonical coordinates. The limit of infinite density and the thermodynamic limit are then 

briefly discussed. 
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I. INTRODUCTION 

^ , The commonly accepted definition of quantum chaos is based on universal statistical properties of suitably defined 
' fluctuations in the energy spectrum. For instance, a confined quantum system with a finite number of degrees of 
OO ^ freedom, e.g., a particle in a billiard, is said to be chaotic whenever the nearest neighbor level spacing (NNLS) 
• distribution, with the spacings normalized to their local average value, is well approximated by the corresponding 
' Wigner distribution for random matrices |Q| . 

. This definition of quantum chaos was conjectured |^ to satisfy a correspondence principle with the dynamical 
Qv^ ' chaos of the associated classical system, i.e., the exponential sensitivity of the classical trajectories to a variation of 
"p*^ , the initial conditions. In fact, a large collection of numerical examples |^ shows that, whenever the classical system 
has positive (non positive) maximum Lyapunov exponent, the corresponding quantum system has Wigner-like (non 
^ \ Wigner-like) NNLS distribution. The recent result gives a theoretical support to this conjecture. 
I ' In the case of systems made of TV identical particles nonlinearly interacting, quantum chaos, in the sense stated 
' O \ above, seems a general rule [^-0] . For these systems, we suggest to look for a correspondence of the quantum chaos with 
^ the dynamical chaos of associated c-number canonical coordinates well distinguished from the classical ones. These 
^ c-number canonical coordinates are appropriate combinations of time-dependent mean fields which approximate the 
dynamics of the A^-body symmetrized or antisymmetrized wavefunction The dynamical equations of the mean 
fields are, in general, nonlinear and allow the presence of dynamical chaos as in classical mechanics |9j. A second 
quantization transforms the mean fields into field operators and restores the exact dynamics of the quantum A^-body 
; ^ ■ system. This second quantization can be made formally identical to the quantization of classical canonical coordinates 
^ and, by analogy, we can expect a correspondence between quantum chaos of A'^-body systems and dynamical chaos of 
their mean field approximations. 

In this paper, we numerically investigate a system of N bosons with both periodic and Dirichlet boundary conditions. 
In Section II, we analyze the exact system from the point of view of quantum chaos and find that a Wigner-like NNLS 
distribution is a general feature in presence of nonlinear interaction already for very low values of N. In Section III, 
we study time-dependent mean field approximations of the same system from the point of view of the exponential 
sensitivity to the initial conditions. We find that the mean fields show dynamical chaos in correspondence to the 
quantum chaos of the exact system for all the values of N considered. We naturally assume that this correspondence 
continues to hold as N increases. In Section IV we quantize the mean field and in in the last Section we discuss the 
limit of infinite density and the thermodynamic limit. 
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II. THE MODEL AND ITS PROPERTIES 



Let us consider a system of N spinless bosons of charge q moving in a one-dimensional lattice with L sites and 
described by the Hamiltonian 
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where index correspondence j ± L = j is assumed. The operator creates a boson in the site j and ctj, and 
7j are the site, hopping, and interaction energies, respectively. Periodic and Dirichlet boundary conditions will be 
considered. In the first case, the system represents a ring threaded by a line of magnetic flux (j) and the phase factors 
are = 2Tr(l)/(j)oL, where 0o = hc/q is the flux quantum (in Gauss electromagnetic system). In the second case, 
the sites lie on a segment and we put /3l = and 9 = 0. The system has wide interest. Its time-dependent 
mean-field approximation has applications to molecular dynamics and nonlinear optics and to electron transport 
in hcterostructures pi] ]. 

All properties of the system can be evaluated by knowing eigenvalues and eigenvectors of H. It is simple to work 
in the space spanned by the Fock states \n\ - ■ ■ n\), where is the number of bosons in the site j and X]j'=i = ^■ 
The index i runs from 1 to the Fock dimension 

N\ (L-l)\ ^ > 

obtained by counting all the possible arrangements of the N identical bosons into the L sites. The D-dimensional 
matrix H representing the Hamiltonian (|^) in the Fock basis has matrix elements 



Hm = (nj • • • n'l\H\n\ • • • n^) = ^ [a^-nj -|- -t^n) {n) - l)] 4i 



1j p 



(3) 



where 



^ , f 1 if n) = n) + 1, nj^+i = nj-^^ - 1, and n\ = nj for / ^ j, j + 1, 

' [0 otherwise. ^ ' 

The eigenvalues Ei and eigenvectors \Ei) of the Hermitian matrix can be numerically evaluated with standard 
methods | [T^ . A bound to the maximum dimension D that can be studied is essentially fixed only by the computer 
memory necessary to storage the full sparse matrix (^). 

For a general, asymmetric system, the NNLS distribution is evaluated from all the eigenvalues Ei . The normalized 
spacings between nearest neighbor levels, whose distribution P(s) is of interest, are taken as 



{E,^r-Ei) I^E^^{i), (5) 



where 



Ai;avW= {Er+k+1 - E,+k) (6) 

with 1 <^ A^av ^ i^- In the case of Dirichlet boundary conditions, the Hamiltonian matrix is real symmetric and 
correspondence to the Gaussian orthogonal ensemble (GOE) of random matrices may be expected. In the case of 
periodic boundary conditions, the Hamiltonian matrix is complex Hermitian and, in general, correspondence to the 
Gaussian unitary ensemble (GUE) may be expected. However, if the external potential, represented by the site 
energies aj , is symmetric under a reflection with respect to some diameter of the ring, the system has an anti- unitary 
symmetry and GOE behavior is restored p^ ]. 

In presence of geometrical symmetries, the level statistics analysis probing the phenomenon of level repulsion is 
meaningful only once the trivial crossings of eigenvalues belonging to different symmetry classes are avoided 
This amounts to analyze the NNLS distribution separately inside each one of the diagonal blocks which compose the 
matrix _ff in a proper basis. For example, a uniform system, i.e., a system with aj, Pj, and "fj independent of the 
site-index j, is invariant under rotation 

n -.j^j + i (7) 
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in the case of periodic boundary conditions, and space- inversion 

V -.j^L + l-j 



(8) 



in the case of Dirichlet boundary conditions. The eigenvalues of H must be divided into L classes, corresponding 
to the eigenvalues exp(i27ri//L), with i/ = of TZ, in the first case, and into two classes, corresponding to 

the eigenvalues ±1 of V, in the second one. This is accomplished by evaluating the operator H in the basis of 
the degenerate eigenvectors oi TZ or V in which the corresponding matrix H is block-diagonal. Each block is then 
independently diagonalised to find the eigenvalues of H within the corresponding symmetry class. 

The eigenvectors li^q) of 7?,, q being the degeneracy index, are obtained by numerically solving the eigenvalue problem 

D 

J2 (J^M - K5m) (nl • • • nl\vq) = (9) 
1=1 

where 

= ("f-4|KK.-nl) . { J - (10, 

The block of H corresponding to the eigenvalue — exp(i27ri^/L) has matrix elements 

D 

= ('"'jl^^ ■ ■■ni)Hk^{n{ ■ ■■nWvp). (11) 

fc,i=l 

A similar, general procedure could be applied also to V . However, the eigenvectors of V are, by inspection, single 
Fock states or symmetric and antisymmetric combinations of couples of Fock states. The even and odd blocks of iJ, 
corresponding to the eigenvalues ±1 of V , have matrix elements which are straightforward combinations of those in 

(i)-. 

Figure 1 shows the NNLS distribution obtained for a uniform system with periodic and Dirichlet boundary condi- 
tions. The distributions of the normalized spacings (H) evaluated for each symmetry class of eigenvalues have been 
summed up for increasing the statistical confidence. The agreement of the calculated NNLS distribution with the 
Wigner surmise for the GOE distribution Pgoe(s) = exp(— •7rs^/4), also shown in the same figure, is statistically 

reliable. 

The importance of performing the level statistics analysis within the appropriate symmetry classes is evidenced in 
Fig. 2, where the NNLS distribution for the same uniform systems of Fig. 1 is evaluated from the total spectrum 
of H. In the case of Dirichlet boundary conditions, the mixing of the even- and odd-parity eigenvalues generates a 
two-peaks distribution which behaves like exp(— s) at large s. In the case of periodic boundary conditions, due to 
the higher number of symmetry classes and their statistical independence, we observe a distribution which mimics 
the Poisson distribution Pp{s) = exp(— s) even at small s |l4j. This fact and the observation that similar results are 
obtained when the symmetries are only approximate |^], will be relevant in the following. 

The above mentioned geometrical symmetries can be explicitly broken by choosing an appropriate j-dependence in 
the parameters a^, /3j, and jj. Figure 3 shows the NNLS distribution obtained with periodic boundary conditions 
when the site energies have the form aj — 2{N — 1)77 ^j, where are arbitrary positive numbers with X]f=i ~ ^^ 
and all the other parameters are as in the uniform case. The resemblance of the calculated distribution with the GUE 
Wigner surmise, -Pgue(s) = (32s^/7r^) exp(— 4s^/7r), expected on the base of the complex Hermitian nature of H, is 
statistically reliable. However, the GOE behavior is restored, as shown in Fig. 4, by choosing symmetric with 
respect to an arbitrary jo- Figures 5 and 6 show the NNLS distribution obtained when the hopping and interaction 
energies have the form indicated in the captions and the other parameters are as in the uniform case. The calculated 
distribution is always close to the GOE one for both periodic and Dirichlet boundary conditions. 

Results similar to those discussed above are obtained for other choices of the parameters, N, L, aj, f3j, jj, and (j) 
i.e., the system (|l|) is generally characterized by a Wigner-like NNLS distribution. Integrability points are the only 
exception to this rule. A first, trivial, point of integrability of the system is the noninteracting case obtained for 
jj — > with eigenvalues given by Ei = J2^=i ^i*^}' being the L eigenvalues obtained for = 1. A second point is 
approached for (3j — > with eigenvalues given by Ei = X^j^i '^j'^} + ^J^K'^i ~ ^ third point of integrability is 
obtained by taking the continuum limit in which Eq. (|l|) becomes the second quantization version of the Hamiltonian 
of A^ bosons interacting via a (5-function potential 
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11=1 



N 

2m dx'i ^ 

" n'=n+l 



c 



6 {xji Xji' ) 



(12) 



with < 
a, = 2/3, 



< 



£. This system is solvable by the Bethe ansatz and can be obtained from (|l|) by putting 
h'^/{mAx'^) and = C/Ax with Aa; = £/L and letting L — > oo. When approaching an integrability 
point, the NNLS distribution transforms into a non Wigner-like distribution whose shape strongly depends on the 
values of the system parameters. 

The results of the level statistics analysis obtained for the boson system (|l]) are in agreement with those found in 
for fermion systems. Quantum chaos, in the sense stated in the Introduction, is a generic feature of systems with 
many particles nonlinearly interacting. The only apparent exception to this rule is a result of [Q| in which a system 
of spinless electrons moving in a ring similar to ours is shown to have Poisson-like NNLS distribution whenever the 
interaction is limited to a small region of the ring. However, studying the same fermion system we found that the 
case considered in Q has an approximate symmetry |^6| . When this symmetry is taken into account or is removed by 
changing the relevant parameters, e.g., adding a random energy to the sites, the NNLS distribution turns to a GOE 
distribution as in the example of Fig. 6. 



III. MEAN FIELD APPROXIMATION 



Two substantially equivalent time-dependent mean field approximations of the Hamiltonian (|^) are obtained by 
choosing the N particles to be described by the normalized boson condensate 



\Z^{t))^iNl)-'/^[al{t)f |0) 



or the normalized coherent state 



\ZM{t)) ^exp y/Nal{t)~y/Na^{t) |0). 
Here, |0) is the vacuum state and a|(t) creates a boson in the single-particle state 

L L 

\z{t)) = alim =Y,z,{t)a]\0) =Y,z,{t)\j). 



(13) 



(14) 



(15) 



The normalization condition of this state, (z(t)|z(t)) — X]^=i — 1; fixes the expectation number of particles in 

the boson condensate (^3|) or the coherent state (^4|) to since in both cases we have (Zjv(t)|a^aj |ZAr(t)) = N\zj{t)\'^. 
The time evolution of the complex amplitudes Zj{t) — {j\z{t)) is obtained from a variational principle for the Dirac 
action M 



dt {ZNmh--H\ZNit)). 

at 



For a general Hamiltonian H = T + V , sum of a single-particle term T and a two-particle term V, 



H = 'y^Tkndj.an + - ^ Vkk'nn' 



kn 



kk'nn' 



we have 



{ZNit)\lh^ - H\ZNit)) = lhJ2NZkit)-^Mt) - ^T^fcn^ Zfc(t)z„(t) 



kn 



I J2 Vkk'nn'N{N ^l)zkit)zk'{t)Zn{t)Zn'{t). 



kk'nn' 



(16) 



(17) 



(18) 



In the above expression as well as in the following ones of this Section, we use the boson condensate (13). Similar 
expressions hold for the coherent state ( p^ with the substitution (iV — 1) —f N. The action ( p^ ) is stationary with 
respect to a variation of Jj (t) if 
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where 

h,i[z{t)] =. T,i + (iV - 1) ^ V.kin Zk{t)z,,{t) (20) 

kn 

are the matrix elements of the mean field single-particle Hamiltonian /i[z(t)]. In the case of the Hamiltonian (|l]), we 
have 

Tji = 6,1 ~ Pie'" d,i+i - ^.e-'" d,+u, (21) 

Vjkin = 27j djk Ski Sin (22) 
and the mean field single-particle Hamiltonian has matrix elements 

hji[zit)] = a.Sji - Pie'^Sj^i+i - Pje~'%j+i.i + 2{N - l)^j\zj{t)fSji (23) 

with j,l = 1, . . . L. 

The system of nonlinear Schrodinger equations (^9|) with hji[z{t)] given by ( p3|) must be numerically solved starting 
from initial values Zj{0) with the normalization condition X^j^i = 1- The choice of the numerical algorithm is 

critically related to the existence of conservation laws. By inspection, the system ( p^ has two conserved quantities, 
the single-particle probability 

L 

\\z{t)f^J2\',{t)f (24) 

and the single-particle energy 

L 

£[z,z] = ^{a,|z,(t)|2 - [/3,_ie^^z,_i(i) +/3je-*«z,+i(t)] z,{t) + {N - l)7,|z,(i)|^}. (25) 

The conservation law (^^ suggests the use a finite-difference Crank-Nicholson scheme 

J2 Ui+'^'^A4t + ^t)]\ z,(t + Ai) =^ ("(5,, -i|J;i,;[z(<)]] zi{t) (26) 
1=1 ^ ^ /=i \ ^ 

which is simple to handle numerically due to the tridiagonal nature of the matrix h . The above scheme would be 
correct C'(At^) if the matrix h were time independent. However, this is not our case and we have to approximate the 



matrix elements in the l.h.s. of (26). The simple approximation hji[z{t + At)] ~ hji[z(t)] has catastrophic effects for 
the conservation law ( p5[ ) unless very small values of At are chosen. An improvement is obtained with an iterative 
procedure. Let us suppose Zj{t) known and try 

Zj{t + At) ^ lim z'f\t + At). (27) 

The n = term is defined by solving 

1=1 ^ ^ 1=1 ^ ^ 

and the n > 1 terms are chosen as solutions of 

L 



E ( + ^^/^./[^^""'Hi + A<)]] z\''\t + At) = E {^^^ - ^^^flW)\ ) (29) 



1=1 ^ ' 1=1 
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This iterative scheme converges in very few steps and ensures an excellent conservation of both quantities (24) and 
The system (|l|) can show local exponential instability. The corresponding maximum Lyapunov exponent A defined 



by 



A=limA(t), A(.).lln^ 
t^oo t ||dz(0)| 



in-5z,{t) = ^/i,,[z(t)]5z,(i) +P[^,(0,<5^j(i),<5^jW] , 



1=1 



where 



(30) 



measures the exponential separation between two states \z{t)) and \z{t)) + e\5z{t)) infinitesimally close (e 0). The 
site projections Szj{t) — {j\5z{t)) satisfy the set of equations obtained by Hnearizing ( |I^ ) 



(31) 



P[z,it),Sz,{t),Sz,it)] - 2(iV~ l)7,z,(t) [z,it)Sz,it) + z,{t)Sz,{t)] 



(32) 



and can be evaluated by numerically integrating ( ^ ) simultaneously with (|T9|). The initial values Szj{0) can not be 
chosen arbitrarily. Indeed, the state \z{t)) + e\6z{t)) must be normalized up to terms 0{e^). Since ||2(t)|| = 1, we 
must have Ke{z{t)\6z{t)) = at any time. However, by using (|l^) and ( ^l| ) we have 



d 



L r 

Re{z{t)\6z{t)) = Re^ 



6zj{t)-^^zjit) + z,{t)^jz,it) 



dt 



= 



(33) 



and, therefore, it is sufficient to have Ke(z(Q)\6z(0)) — 0. 

For solving (|3l] ) simultaneously with ([l^) we again adopt an iterative modification of the Crank-Nicholson scheme 
in which P is considered as a driving term. Let us suppose Zj (t) and dzj (t) known and try 



Szj(t + At)= lim (5z]"\i + At) 



(34) 



The n = term is defined by solving 



1=1 



E ( ^si + ^^h,i[z(t)] ) dzf\t + At)=J2( 6,1 - t—h.Mt)] ) 5zi{t) 



1=1 



.At. 



2h 



-i^{P [^J it) ' (*) ' i^J it) ' it) ' it)] 



and the n > 1 terms are chosen as solution of 
^ ' .At. 



1=1 



E ( + ^^h,,[z("-'Ht)] ) 5z\-\t + At)=Y,[ h ^^h,,[zit)] ) 6zi{t) 



1=1 



.At. 



2h 



.At 



(f + At),Sz^''~'\t + At), 5zf-'-\t + At)] + P[z,{t),6z,{t),6zj{t)]^ 



(35) 



(36) 



{t + At) and (5zj""^'(i + At) 



Note that (pq) must be solved after (M) and (36) have been solved at the iteration n — 1 in order to know both 



The quantity ||(5z(t)|| can grow exponentially and, therefore, Szj[t) must be periodically scaled in order to avoid 
numerical overflows IQ. The scaling factors are stored for computing the Lyapunov exponent (^0|). When the 
system ( |l9|) is chaotic, the computer round-off errors inevitably make the numerical solutions obtained with different 
integration steps At different after a sufficiently long time. Therefore, the comparison of solutions relative to different 
steps is not a good check that the algorithm correctly works, unless time-averaged quantities, e.g., the Lyapunov 
exponent (|30|), are compared [|l9|. On the other hand, a check based on the conservation of the quantities (^4|) and 
( pq ) is meaningful and can be used to fix the size of the integration step in relation to a chosen accuracy. With the 
modified Crank-Nicholson scheme described above, for At ^ 10~^h/r] after 10^ iterations we have relative errors in 
(p3) and (|25|) which are smaller than 10^^ and 10^^, respectively. 
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Figures 7 and 8 show the behavior of A(t) with periodic and Dirichlet boundary conditions, respectively. The curves 
denoted with u, a, (3, and 7, refer to the system in which all the parameters are independent of j as in Fig. 1, Uj is 
random as in Fig. 3, (3j is as in Fig. 5, and 7^ is localized as in Fig. 6, respectively. After an initial transient, not 
shown in Figs. 7 and 8, A(t) approximately stabilizes around a positive value which we take as the corresponding 
maximum Lyapunov exponent A. Note that A~^ < lO^h/r] is much smaller than the maximum simulation time lO^h/r]. 
The value of A is independent of changes in the initial conditions 6zj{Q). It is also independent of changes in the 
initial conditions 2:^(0) provided that the conserved energy £[z,'z] is not changed. The maximum Lyapunov exponent 
exceptionally vanishes when the initial state 1^(0)) is coincident or very close to one of the stationary states of the 
system (|l9|) which are defined by 

\zE{t))^e-i^'\zEm, \\zE{t)r = ^- (37) 

The comparison of Figs. 1-6 with Figs. 7-8 suggests a correspondence between quantum chaos of a system of 
interacting particles and dynamical chaos of its mean field approximations. Whenever the exact system (n) shows 
Wigner-like NNLS distribution the corresponding mean field system ([l9|), or that for the coherent state (^), has a 
positive maximum Lyapunov exponent. 

The simultaneous presence of chaotic behavior in the exact and mean-field systems is obtained also for values of the 
parameters iV, L, aj, Pj, 7j, and (f) different from those reported. In the three cases in which the exact system has 
been shown to be integrable, namely jj — > 0, Pj — > 0, and the continuum limit, the corresponding mean field system 
is integrable and has A = 0. For 7^ — > 0, the system (^9|) becomes linear and A{t) = 0. For (3j — > 0, we have 



and, therefore, 



^ |z, (t) p = z, (t) ^-z, it) + z, it)j^ zj (t) = (38) 



z^it) = z^.(0)e-^(«^+2(^-i)^.l-.(o)l')*A. (39) 

The corresponding variation 

-.,(0).i2(iV-lh,|i.,(0).,,0, + .,(OM.,(0)le-<".«"-'>-*.m|-../« ,40, 

shows that |l(5z(i)|| is 0{t) and therefore A = 0. Finally, in the continuum limit the system ( p^ ) becomes the well 
known nonlinear Schrodinger equation, 

in^^z{x, t) = -^-^z{x, t) + C{N - l)\z{x, t)\^z{x, t), (41) 
solvable via spectral transform [pO[. 



IV. SECOND QUANTIZATION OF THE MEAN FIELD 

The correspondence between quantum chaos of an exact A^-body system and dynamical chaos of the associated mean 
field approximations parallels the correspondence between quantum chaos of a single-particle system and dynamical 
chaos of the associated classical equations. This parallel can be connected to the fact that the first quantization of 
the classical canonical coordinates has a formal counterpart in the second quantization of the mean field Let us 
see this in detail. The mean field Zj{t) ~ ^/Nzj{t), j = 1, . . . , L, is determined by the dynamical equation 

ihj^Zjit) - ajZ.j{t) - /3,_ie'''Z,_i(i) - pjC-'" Z^+^it) + 27,Z,(i)Z,(t)Z,(i). (42) 

In this Section we will consider the coherent state ( [T^ ) but similar results hold for the boson condensate (|l3| ) with 
the substitution N [N — 1). The quantization rule 

Z,{t)^Z,{t), Z,{t)^Z]{t), (43) 

with 
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0, 



Sjk 



(44) 



transforms the nonhnear Schrodinger equation ( [4 21) into the Heisenberg equation for the field operator Zj{t). Indeed, 

in the representation of the site- localized states (j)'- = 6jk, we have Zj{t) — J2k=i 4''j'^k{t) — o,j{t) whose Heisenberg 
equation of motion is 



(45) 



The second quantization (^^ can be made formally identical to the first quantization of classical canonical coordi- 
nates with the standard transformation 



1 h 



Qj (0 = \/ o ( (t) + z, it) ) , p,- (i) = - y - ( Z, (t) - Z, it) 



By using the total energy of the system 

n[Q,P] = N£[z{Q,P)XQ,P)] 



2fi 



^,eQ,-i{t)+iP,-i{t) 



2n 



2h 



,{t)~lpj(t) 



2% 



Q,{tf + P,{tr ^^ 
2h 



(46) 



(47) 



the nonlinear Schrodinger equation (^) can be rewritten in the Hamilton formalism 



dPj{t) 



■H[Q,P], -P,{t) 



dQj{t) 



■H[Q,P] 



(48) 



and the quantization (|44| ) is equivalent to the introduction of Hermitian operators Qj (t) and Pj {t) with commutation 
relations 



Q,(t),Qfe(<) = Pj{t).Pk{t) 



0: 



l{t),Pk{t) 



ihSjk- 



(49) 



V. THE LIMIT N ^oo 



We first analyze the limit of infinite density, in which iV — s- oo with L constant. 

In this limit the quantum theory reduces to a c-number theory and is, therefore, analogous to a classical limit even 
though fi ^ Q [^^ . For our system this is easily seen by retracing the N dependence in the equations of the previous 
section. For large N , the system reduces to a collection of independent nonlinear oscillators whose nonlinearity grows 
like A'^. From the point of view of chaotic properties, the limit N / L cx3 is therefore trivial. An example where the 
same kind of limit gives rise to a non trivial chaotic system can be found in . 

It is interesting to see the emergence of the limiting behavior N / L — > oo by comparing the cumulative density of 
states 

D 

D{E) = Tr d{E - H) = J29{E ~ E,) (50) 

with the approximate expression obtained according to the Weyl rule 

Anf(^) = J dQi---dQL J dPi---dPL 9{E-n[Q,P]) diN-M[Q,P]). (51) 

Note that we have a (5-function constraint on the Q — P phase space which fixes 
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^ Q,(t)2 + P,(t) 



U[Q,P]^Y. ^''> =N (52) 



in agreement with (24). Due to the presence of this constraint, the r.h.s. of ( pl| ) is a {2L— l)-niultiple integral which 
can be evaluated with Monte Carlo integration in the hypercube of side 2\/2^JV centered in the origin. Figures 9-11 
show, in a case with Dirichlet boundary conditions, that when the density p = N / L is increased D{E) is approximated 
by Dm[{E) with increasing precision. Similar behavior is obtained with periodic boundary conditions and/or different 
values of the parameters of the system. The total number of levels given by ( pO[ ) and ( |5l| ) become equal in the limit 
N/L — > oo. Indeed, we have 

D^t = Yiu\^ Dai{{E) = J dxi... dx2L ^ j = _ (53) 

which is the value of the Fock dimension (||) for N ^ L. 

For a single-particle system, the smooth behavior of the cumulative density of states is approximated by the 
corresponding semiclassical expression in the limit of high energies (?i ^ 0). Analogously, in the case of an A'^-body 
system for N/L large we can use Da\f{E) to approximate the smooth behavior of D{E) and evaluate the normalized 
spacings (||) according to 

= {E,+i - E,)^D^i{E) ~ D^i{E,+i) - D^i{E,), (54) 
dE 

as suggested in [|] . 

Finally, let us briefly discuss the thermodynamic limit. Unlike the limit N/L — ^ oo, the system preserves its 
quantum features and the mean fields do not give a complete description. This is also reflected by the behavior of the 
cumulative density of states. For the system discussed here, when N and L ^ oo with N/L = p constant by using 
( p3| ) and (||) we have 

InAnf _ 1 + 

^ (l+p)ln(l + p)-plnp - ^^^> 

The function p{p) is smaller than unity for any finite p and tends to unity for p — oo. Therefore, I?,nf/£) « iP^Cp)-! 
vanishes in the thermodynamic limit. 

On the base of our numerical results and the considerations made in the previous Section, the correspondence 
between quantum chaos of an iV-body system and dynamical chaos of its mean-field approximations can be naturally 
assumed to hold for N ^ oo. This fact is exploited in |2^ where the authors consider the chaotic behavior of the 
same system discussed here when the thermodynamic limit is approached. 
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FIG. 1. NNLS distribution P{s) for the system (|l[) in the uniform case N — 5, L = 9, Oj = 0, f3j — r;, and = r/ 
with 0/00 = 0.3 and periodic boundary conditions (solid histogram) and (j>/4'o ~ and Dirichlet boundary conditions (dashed 
histogram). The solid line is the Wigner surmise for the GOE distribution. 
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FIG. 2. As in Fig. 1 but without separating the eigenvalues into the appropriate symmetry classes. The dashed line is the 
Poisson distribution. 
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FIG. 3. As in Fig. 1 with periodic boundary conditions but choosing aj = 2{N — l)r] ^j, where are arbitrary positive 
numbers with X^^^j = 1- The soUd and dashed hnes are the Wigner surmise for the GOE and GUE distributions, respectively. 
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FIG. 4. As in Fig. 3 but clioosing ^jo-j = Cjo+i with jo arbitrary. 
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FIG. 5. As in Fig. 1 but choosing (32 = Ps = 0.5r; and Pj = ?? for j / 2, 3. 
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FIG. 6. As in Fig. 1 but choosing 7^- = rjSjs. 
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FIG. 7. Maximum Lyapunov exponent for the mean field system with periodic boundary conditions. The curves 
denoted with u, a, f3, and 7, are obtained with the parameters given in Figs. 1, 3, 5, and 6, respectively. The initial conditions 
ZjijS) and SzjijS) are a set of arbitrary complex numbers satisfying ||z(0)|p = 1 and Yle{z{Q)\5z{Q)) = 0. 
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FIG. 8. As in Fig. 7 with Dirichlet boundary conditions. 
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FIG. 9. Cumulative density of states in the case A'' = 10, L = 3, aj, 13 j, and 7j chosen between and 77 with a random 
number generator, and with (j>/(j)o ~ and Dirichlet boundary conditions. The histogram is the exact result and the 
dashed line is the mean field approximation (pi]). 
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FIG. 10. As in Fig. 9 with N = 30. 
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FIG. 11. As in Fig. 9 with N = 90. 
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